#include "GridOperator.cuh"
#include "GridSampler.cuh"
#include "GridPainter.cuh"

template<typename T>
__global__ void cylinderKernel(int nx, int ny, int nz, int z, T dx, T cx, T cy, T cz, T r, T h, T v, T *f)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;
	if(x < nx && y < ny)
	{
		int idx = x + y * nx + z * nx * ny;
		T dcx = (x + T(0.5)) * dx - cx;
		T dcy = (y + T(0.5)) * dx - cy;
		T dcz = (z + T(0.5)) * dx - cz;
		if(dcx * dcx + dcz * dcz <= r * r && abs(dcy) <= h)
			f[idx] = v;
	}
}

template<typename T>
void cylinder(int nx, int ny, int nz, T dx, T cx, T cy, T cz, T r, T h, T v, T *f)
{
	dim3 gdim((nx - 1) / BLOCK_DIM + 1, (ny - 1) / BLOCK_DIM + 1);
	dim3 bdim(BLOCK_DIM, BLOCK_DIM, 1);
	for(int z = 0; z < nz; ++z)
		cylinderKernel<T><<<gdim, bdim>>>(nx, ny, nz, z, dx, cx, cy, cz, r, h, v, f);
	cutilCheckMsg("cylinderKernel() execution failed!");
}

template
void cylinder<float>(int nx, int ny, int nz, float dx, float cx, float cy, float cz, float r, float h, float v, float *f);

template<typename T>
__global__ void sphereKernel(int nx, int ny, int nz, int z, T dx, T cx, T cy, T cz, T r, T v, T *f)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;
	if(x < nx && y < ny)
	{
		int idx = x + y * nx + z * nx * ny;
		T dcx = (x + T(0.5)) * dx - cx;
		T dcy = (y + T(0.5)) * dx - cy;
		T dcz = (z + T(0.5)) * dx - cz;
		if(dcx * dcx + dcy * dcy + dcz * dcz <= r * r)
			f[idx] = v;
	}
}

template<typename T>
void sphere(int nx, int ny, int nz, T dx, T cx, T cy, T cz, T r, T v, T *f)
{
	dim3 gdim((nx - 1) / BLOCK_DIM + 1, (ny - 1) / BLOCK_DIM + 1);
	dim3 bdim(BLOCK_DIM, BLOCK_DIM, 1);
	for(int z = 0; z < nz; ++z)
		sphereKernel<T><<<gdim, bdim>>>(nx, ny, nz, z, dx, cx, cy, cz, r, v, f);
	cutilCheckMsg("sphereKernel() execution failed!");
}

template
void sphere<float>(int nx, int ny, int nz, float dx, float cx, float cy, float cz, float r, float v, float *f);

template<typename T>
__global__ void rectangleKernel(int nx, int ny, int nz, int z, T dx, T cx, T cy, T cz, T w, T h, T d, T v, T *f)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;
	if(x < nx && y < ny)
	{
		int idx = x + y * nx + z * nx * ny;
		T dcx = abs((x + T(0.5)) * dx - cx);
		T dcy = abs((y + T(0.5)) * dx - cy);
		T dcz = abs((z + T(0.5)) * dx - cz);
		if(dcx * T(2) <= w && dcy * T(2) <= h && dcz * T(2) <= d)
			f[idx] = v;
	}
}

template<typename T>
void rectangle(int nx, int ny, int nz, T dx, T cx, T cy, T cz, T w, T h, T d, T v, T *f)
{
	dim3 gdim((nx - 1) / BLOCK_DIM + 1, (ny - 1) / BLOCK_DIM + 1);
	dim3 bdim(BLOCK_DIM, BLOCK_DIM, 1);
	for(int z = 0; z < nz; ++z)
		rectangleKernel<T><<<gdim, bdim>>>(nx, ny, nz, z, dx, cx, cy, cz, w, h, d, v, f);
	cutilCheckMsg("rectangleKernel() execution failed!");
}

template
void rectangle<float>(int nx, int ny, int nz, float dx, float cx, float cy, float cz, float w, float h, float d, float v, float *f);
